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NUMERICAL CALCULATIONS OF VELOCITY AND PRESSURE 
DISTRIBUTION AROUND OSCILLATING AIRFOILS 


Theodore Bratanow,* Akin Ecer,** 
and Michael Kobisket 
University of Wisconsin-Milwaukee 


SUMMARY 


An analytical procedure based on the Navier-Stokes equations was devel- 
oped for analyzing and representing properties of unsteady viscous flow 
around oscillating obstacles. A variational formulation of the vorticity 
transport equation was discretized in finite element form and integrated nu- 
merically. At each time step of the numerical integration, the velocity 
field around the obstacle was determined for the instantaneous vorticity 
distribution from the finite element solution of Poisson's equation. The 
time-dependent boundary conditions around the oscillating obstacle were in- 
troduced as external constraints, using the Lagrangian Multiplier Technique, 
at each time step of the numerical integration. The procedure was then ap- 
plied for determining pressures around obstacles oscillating in unsteady 
flow. The obtained results for a cylinder and an airfoil were illustrated 
in the form of streamlines and vorticity and pressure distributions. 


INTRODUCTION 


The objectives of the investigation were to develop a method for analy- 
sis of general, unsteady, viscous and incompressible flow around an oscil- 
lating obstacle of arbitrary shape and to determine the velocity and pres- 
sure distribution around this obstacle. Using the finite element method, a 
numerical procedure was developed and tested by investigating first the un- 
steady flow past a stationary circular cylinder. This method was then ap- 
plied to analyze the flow around stationary and pitching airfoils. The vel- 
ocity and pressure distributions around these obstacles were determined. 

The method of analysis considers the viscous effects of the boundary layer 
and the overall flow characteristics in a unified mathematical frame. 

Unsteady flow around a circular cylinder has been analyzed by several in- 
vestigators [1,2, 3, 4, 5]. In most cases, however, the numerical procedure 
was developed specifically to suit the geometry of a stationary obstacle. 

The problem of moving obstacles was analyzed only for a circular cylinder os- 
cillating around its longitudinal axis [4] , where the instantaneous geometry 
of the obstacle remains the same. A review of existing studies of unsteady. 
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viscous flow showed discrepancies in the results between analyses and exper- 
iments. In the case of a circular cylinder, the experimental results for 
the surface pressure distribution do not agree with theoretical results [3, 
6,7]. Similar discrepancies are observed for the drag coefficients and the 
points of separation. 


SYMBOLS 
2 2 

A area of a finite element, in. (m ) 

A rectangular matrix defined in Eq. (Bl) 

A general differential operator 

^ rectangular matrix representing the boundary conditions on the obstacle 

rectangular matrix containing matrix b* = 

^ row vector defining the variation of vorticity, pressure and velocit- 
ies over an element 

£ rectangular matrix defining the free stream boundary conditions 

£ row vector defining the variation of stream functions over an element 

£ column vector specifying the stream function values on the free stream 
boundary 

£ column vector specifying the vorticity on the obstacle's boundary 

f f(x,y) - arbitrary function constrained by boundary conditions 

F F(x,y,A^3 - arbitrary function with Lagrange multipliers as defined 
in Eq. (27) 

£ inverse of matrix £* 

h column vector defined in Eqs. (49), (57), and (65) 

h* column vector containing vector h_, li* = (li £) 

£ identity matrix 

£ rectangular matrix defined in Eq. (B6) 

k 
K 


2 


constant members in Eqs. (B5) and (B9) 

row vectors representing the stream function's partial derivatives 
with respect to the area coordinates 



-M. M(UjV,(jJ) - nonlinear term in the vorticity transport equation, 

M = u + V 

9x 9y 

£ null matrix 

2 2 

p local static pressure, Ibf/in. (N/m ) 

£ column vector representing the pressure values at the nodal points of 
an element, lb£/in.^ (N/m^) 

£ column vector representing the pressure values at the nodes of the en- 
tire grid, Ibf/in.^ (N/ra^) 

£* column_vector containing vector £ and the Lagrange multipliers, 

r = ce 

£ symmetric matrix defined in Eq. (65) 

£* symmetric matrix containing matrix £ and the free stream boundary 
conditions 


E 

Q 


column vector resulting from the multiplication of column h* by a 
symmetric inverse 


Q(x,y), Q = 2 


9u 9v 
9y 9x 


9u 9v 
9x 9y 


£ column vector in Eq. (40) 

£ rectangular matrix resulting from the multiplication of matrix £* by 
a symmetric inverse 


£ symmetric matrix defined by Eq. ( 49 ) 

£* symmetric matrix containing matrix £ and the free stream boundary 
conditions 


t time, sec. 

£ square transformation matrix defined in Eq. (B7) 

velocity components, in. /sec. (m/sec.) 

£,£ velocity vectors, in. /sec. (m/sec.) 

3 3 

V volume of a fluid element, in. (m ) 
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w 

w* 

x,y, 

z 

X,Y 

V 

v2 

V 

r 

X 

V 
P 

i 

<p 

i 

i 

0 
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symmetric matrix defined by Eq. (57) 

symmetric matrix containing matrix W and the free stream boundary con- 
ditions 

coordinate directions 

body forces in the x and y directions respectively, Ibf (N) 

9 9 

del operator, ^ 2. 

Laplace operator, j + — j 

3x^ 3y^ 

tangential and radial coordinate directions respectively 
column vector representing the Lagrange multiplier 

column vector representing the second derivatives of the stream func- 
tion at the nodal points of an element 

length of the perimeter of area A, in. (m) 

Lagrange multiplier 

2 2 

kinematic viscosity coefficient, in. /sec. (m /sec.) 

3 3 

density of the fluid, Ibm/in. (kg/m ) 

constant members of the matrix in Eq. (B6) 

i^(x,y) - boundary conditions of ah arbitrary function 

general quadratic functional for variational formulations 

2 2 

’i’(x,y) - stream function defined in Eq. (6), in. /sec. (m /sec.) 

column vector representing stream function and velocity values at the 
nodal points of an element 

column vector representing the stream function and velocity values at 
the nodes of the entire grid 

^olumn_vector containing vector and the Lagrange multipliers, 
i* = ci ^3) 

angle between the horizontal line and the line normal to the surface 
of the obstacle 



w vorticity at a point, sec. 

w vorticity vector, sec. ^ 

0 

u column vector representing the time derivative of vorticity at the 

_2 

nodes of the entire grid, sec. 

• • 

w* column vector containing vector w and the Lagrange multipliers, 

(b* = (w 6 ) 

— ^ cw^ 

C area coordinate, in. (m) 

Subscripts 
b obstacle 

c free stream 

i i-th boundary condition 

m m-th node 

n n-th node 

o at time t 

o 

p pressure 

s stream function 

V total volume 

w vorticity 

x,y derivative with respect to x or y respectively 

5 tangential direction 

ri normal direction 

Superscripts 

variable powers 


t 


transpose 



MATHEMATICAL ANALYSIS 


The mathematical analysis of the two-dimensional unsteady viscous flow 
around an obstacle of arbitrary shape will be described. The variational 
formulation of the problem is based on the most general form of the Navier- 
Stokes equations. The finite element method was applied for the discrete 
representation of the governing equations. The equations from which veloc- 
ity and pressure distributions were calculated, were derived using the prin- 
ciple of conservation of mass and the balance of momentum in terms of stream 
functions and vorticities. The boundary conditions applied for the solution 
of these equations will also be described. 


Equations of Motion 

The governing differential equations of motion for the analysis of two- 
dimensional unsteady flow of incompressible viscous fluids can be written 
as follows [8] : 


3u 

8t 


+ 



9v 

9t 


+ 



= i X - +v 

2 2 
a u . a u 

3y p p 8x 

2 2 
L 9x ay^ J 

^ ly - 

2 2 ~ 
a V 9 V 

ay p p ay 

2 ' 2 
L 9y^J 


CD 


C2) 


9u ^ 
3x 9y 


0 


(3) 


Eqs. (1) and (2) are the Navier-Stokes equations and Eq. [3) is the continu- 
ity condition for incompressible fluids. 


Stream function and vorticity formulation - The vorticity and pressure 
distributions can be obtained from a simultaneous solution of Eqs. (1), (2) 
and (3). Although a direct solution of these equations is possible, the 
equation of continuity can be included in the formulation in a more natural 
manner using the concept of stream functions and vorticities. The stream 
functions can be defined in such a way that they satisfy the equation of 
continuity. The system of Eqs. (1), (2) and (3) can then be reduced to two 
equations in terms of vorticities and stream functions. The vorticity is 
defined as [8] : 

(jj = V X u (4) 
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In the case of two-dimensional flow only one component of the vorticity 
exists and it is perpendicular to the plane of the flow: 


8v 3u 
3x 8y 


(5) 


The stream functions are defined for the two-dimensional case as follows: 


_ 

3y 



( 6 ) 


Substituting Eq. (6) into Eq. (3), the continuity equation is automatically 
satisfied. From Eq. (5) the vorticity can be written in terms of the stream 
functions as follows: 

03 = - (7) 


Vorticity transport equation - Vorticity is a measure of the angular mo- 
mentum of the fluid particles in motion. Eqs. (1), (2) and (3) can be com- 
bined and using Eqs. (4), (5) and (6), the vorticity transport equation or 
the Helmholtz equation is obtained as follows [8] : 


3o3 3o) 

3t ^ 3x 


3o) _ 
3y “ 


V 




C8) 


According to this equation, the time derivative of the vorticity in space, 
which consists of local and convective terms, is equal to the rate of dis- 
sipation of the viscous terms. 

Using Eq. (7), Eq. (8) can be converted into another equation with a 
single variable in terms of the stream functions as follows: 


3\|3 3V^3|3 3\|3 37^\j> 

3t * 3y 3x 3x 3y 




C9) 


The left-hand side of the equation again contains local and convective in- 
ertia terms, while the right-hand side describes the viscous dissipation. 
When the fluid motion is slow, the viscous forces are considerably greater 
than the convective inertia terms. Neglecting the convective inertia terms, 
the equations of motion become linear with respect to velocity and describe 
the creeping motion of a viscous fluid. 


At high Reynolds numbers the inertia forces are more important in com- 
parison with the viscous forces. In this case, the viscous forces can not 
be neglected because, although they are small in magnitude, they change the 
basic characteristics of the equations. If the viscous terms are neglected, 
the complete set of boundary conditions is not necessary for a unique solu- 
tion of the problem. This corresponds to the Euler's equation for an ideal 
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fluid. A unique solution of Euler's equation is obtained by using a bound- 
ary condition specifying zero normal velocity at the boundary. However, the 
condition for zero tangential velocity is not satisfied. Therefore, the 
boundary layer is still important at higher Reynolds numbers although the 
changes from the ideal fluid conditions occur only in a small region around 
the obstacle. In the conducted analysis, the Navier-Stokes equations were 
kept in the complete form for the cases of different Reynolds numbers and 
the boundary layer effects were considered. 


Equation for pressures - The pressures can be calculated from another 
equation, which can be derived again from the Navier-Stokes equations. 

Using Eqs. (1), (2) and (3), the pressure distribution is expressed in terms 
of the velocity distribution as follows [9] : 


-p<i 

3x^ 3y^ 

where 

_ ^ ^ ^ ^ 

^ “ 8y 3x “ 3x 3y 


( 10 ) 


( 11 ) 


According to Eq. (10), the instantaneous pressure distribution can be de- 
termined for a particular velocity distribution. On the other hand, pres- 
sures are not required for the determination of the velocity distribution. 


Boundary Conditions 


The boundary conditions for stream functions, vorticities and pressures 
around the obstacle and the free stream boundary conditions can be des- 
cribed as follows: 


Boundary conditions for stream functions - For the solution of Poisson's 
equation, Eq. (7), two types of boundary conditions were considered: 
Dirichlet type boundary conditions, specifying the stream functions, and 
Neumann type boundary conditions, specifying the normal slope of the stream 
functions. These two types of conditions can be summarized as follows [11]: 

^ = prescribed for free stream flow on the outer rectangular boundary 
= zero on the outer rectangular boundary 


" 5 ^ = free stream velocity on the outer rectangular boundary 
dy 

^ = constant on the obstacle boundary 
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^ = zero on the obstacle boundary (see figure 4) 

os 

Boundary conditions for vorticity - The conditions for the free stream 
boundary can be defined according to Kelvin's theorem. A closed boundary 
in the free stream will have zero vorticity along the boundary line. On the 
other hand, the vorticity at the obstacle boundary can be defined in terms 
of the normal derivative of the tangential velocity only. 

0 ) = zero on the outer rectangular boundary 
9ii 

0 ) = ^ on the obstacle boundary 

9n 

Boundary conditions for pressures - The boundary conditions for pressures 
by specifying the static pressure at the free stream boundary 
the fluid pressure around the solid boundary. 

the outer rectangular boundary 
the obstacle boundary 


can be defined 
and specifying 

p = zero on 

7 ^ = zero on 
9n 


Variational Formulation of 
Navier-Stokes Equations 


A variational method was applied for the solution of Eqs. (7), (8) and 
(10). The solution of the following partial differential equation 

A ip = oi ( 12 ) 


was sought, under a set of boundary conditions, where A represents the dif- 
ferential operator. The solution of Eq. (12) can be calculated by minimiz- 
ing the quadratic functional 


- 2 <I> = 



(13) 


with respect to solution functions that satisfy the boundary conditions. 


Stream function equation - The variational formulation of the Poisson's 
equation is the well-known Dirichlet's problem [10]. As applied to the 
Poisson equation, Eq. (7), Eq. (13) becomes 


2 f = 



(14) 
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The first term of the integral in Eq, .(14) can be rewritten using Green's 
theorem as 



ip dA = - 





dip 

3n 


dr 


(15) 


The last integral in Eq. (15) vanishes because of the applied boundary condi 
tions. One of the boundary conditions for Eq. (7) states that the normal 
derivative 9i|^/8x is zero on the rectangular boundary. Thus, the value of 
the integral is identically zero on the vertical sides of the rectangle. On 
the upper and lower sides, the stream functions have constant magnitudes but 
alternating signs. Thus, their sum is equal to zero. After the application 
of Green's theorem, Eq. (7) can be written as 



The solution of Eq. (7) can now be obtained by minimizing the functional $ 
in Eq. (16). The stream functions and the velocity distribution can be ob- 
tained from this formulation for a given set of vorticities. 


Equation for pressures - The variational formulation of the equation for 
pressures, Eq. (10), is again a Dirichlet type of problem. A similar pro- 
cedure is followed for the variational formulation of the Poisson's equa- 
tion for pressures, Eq. (10). The contour integral in Eq. (15) vanishes 
again, since the pressure is identically zero on the outer boundary. The 
resulting variational functional is 



The parameter Q in Eq. (17) is defined for a particular distribution of the 
velocity. 

Vorticity transport equation - For the finite element representation of 
the vorticity transport equation, Eq. (8), a variational formulation of this 
differential equation was developed. However, for this case, an exact var- 
iational formulation does not exist due to the nonlinear terms on the left- 
hand side of Eq. (8). The right-hand side is in the form of Poisson's equa- 
tion and can be written in a variational form corresponding to Dirichlet 's 
problem. The solution of the right-hand side of Eq. (8) can be obtained by 
minimizing a functional given as 
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dA 


(18) 


$ 


1 


V 


2 





The local time derivative of the vorticity term in Eq. (8) can be formu- 
lated in variational form using Hamilton's principle. The discretized form 
of the vorticity transport equation can then be written in matrix form as 

9w 

^=Mw (19) 


where the matrix M is a nonlinear function of velocity and vorticity dis- 
tributions. The solution of Eq. (19) can be obtained by minimizing the 
functional given as 




dA 


( 20 ) 


The matrix M in Eq. (19) can be derived from a perturbation analysis of the 
vorticity transport equation. From a Taylor series expansion of the vel- 
ocities in terms of the vorticities, the velocity distribution at time t 
can be written as 


u = u + 
— —0 1 

1 

1 9w 1 

(Aw) + 

fd\\ 
1 2 


1 -1 

0 ' 

i 3 w J 


(Aw) 2 

2 ! 


+ 


( 21 ) 


Once the instantaneous derivatives of the velocities with respect to vor- 
ticities are determined, the matrix M can be calculated for any desired 
accuracy. From the equation of continuity and the definition of vorticity 

V*u=0 Vxu = to (22) 


the derivatives of velocities with respect to vorticities can be calculated 
as 


V 


3u 

9(j0 


0 


V 


9u 

9co 


I 


(23) 


The solutions of the systems of equations in Eq. (23) are similar to the 
solution of Eq. (7) except for a right-hand side consisting of unit vortic- 
ities at each nodal point of the finite element grid. Considering that the 
velocity distribution is not very sensitive to the incremental changes in 
the vorticity distribution, only the first terms in the Taylor series can be 
employed in the analysis. The variational formulation of the vorticity 
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transport equation can then be written as follows: 



The local time derivatives of the vorticities can be calculated by minimiz- 
ing Eq. (24) . 


Lagrange Multipliers 


The variational formulation of the governing differential equations de- 
rived above does not include the boundary conditions. The time dependent 
boundary conditions can be introduced in a simple fashion which then becomes 
one of the main advantages of the variational method. The Lagrangian mul- 
tiplier technique was used for this purpose. Accordingly, the boundary con- 
ditions are added to the variational formulation as additional constraint 
conditions. The method of solution will now be summarized. 

Considering a variational functional of the form 

f(x,y) (25) 

where the boundary conditions are excluded and that 

ii(x,y) = 0 (26) 

is a set of boundary conditions, the solution of the differential equations 
satisfying the boundary conditions is sought through another functional 

F(x,y,A^) = f(x,y) ^iij^(x,y) (27) 

where is the Lagrange multiplier. The solution can be obtained by mini- 
mizing F(x,y,A^) with respect to x, y and A^, as 





(28) 


In the case of a quadratic functional, a system of linear algebraic equa- 
tions are obtained. The application of the Lagrange multiplier technique 
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produces a constrained minimum of the functional, which satisfies the bound- 
ary conditions. 

New boundary conditions can be included using the Lagrangian multiplier 
technique by defining only the additional constraint equations at each time 
step. In this problem the obstacle boundary conditions change at each step 
of the numerical integration, depending on the position and the motion of 
the obstacle. For a stationary obstacle, the Poisson’s equation is solved 
only for changing terms on the right-hand side of Eq. (7). The boundary 
conditions at the free stream boundary do not change with time. The appli- 
cation of the Lagrangian multiplier technique reduces the computational cost 
of the solution considerably. 


THE NUMERICAL SOLUTION OF THE GOVERNING EQUATIONS 


When the variational formulation of the governing differential equations 
is completed, they can be discretized in finite element form. The result- 
ing discrete systems consist of linear algebraic equations and a matrix dif- 
ferential equation. The details of the discretization and the application 
of the finite element method is summarized in this section^. The solutions 
of the systems of equations are analyzed and the corresponding boundary 
conditions are discussed. 


Discretization of the Equations of Motion 


Following the basic discretization procedure of the finite element anal- 
ysis, the obstacle and the surrounding continuum were represented by tri- 
angular finite elements. The finite element gridwork for the analysis of 
flow around a cylinder is shown in figure 1. In this case the objective 
of the particular triangulation was to represent most accurately, the 
boundary layer effects around the cylinder. However, in general, it is not 
necessary to describe the obstacle boundary with nodal points. For example, 
the gridwork in figure 2 was used for the analysis of the flow around a 
moving airfoil, where the position of the airfoil changed at each time step. 

Stream function equation - The variation of the stream functions over a 
triangular element was approximated by a third-order polynomial , which can 
uniquely be determined from the values of the stream functions and veloci- 
ties at the nodes of the triangles. The stream function at any point of 
the triangle can then be written as 

ij; = C A (29) 

where A is a rectangular matrix describing the approximation functions as 
shown in Eq. (Bl) of Appendix B. ^ is a row vector formed by the area co- 
ordinates of a point and is another vector formed by the nodal veloci- 
ties and the stream functions of the finite element: 
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^ 

^2 ^3 ^ 1^2 


^ 1^3 ^ 2 ^ ^ 2^1 ^ 3^1 ^ 3^2 ^ 1 ^ 2^3 


(30) 


,t 

ij; = 


3ij; 9t|j 

h ^ W "^2 


‘^A>2 3t|>2 91^3 94^3 

3x 9y *^3 9x 9y 


(31) 


The area coordinates ?. in Eq. (30) are described in detail in Appendix A 

[11]. " 

Vorticities, on the other hand, were considered to vary linearly over a 
triangular element and can be defined as follows: 


w = 


'l «2 %] 


w. 






= B CO 


(32) 


Eqs. (29) and (32) can be substituted in Eq. (16) to obtain a discrete var- 
iational functional as follows: 


* = i 


'k 


,t 

~ t t t t 

ACC A + A C C A 

dA - / 

A^cS 


X— X— — _ 


- 


CO dA (33) 


In this case the integration domain A indicates an integration over the area 
of the triangular element and a summation of the terms at each node from all 
neighboring elements. 

The solution of Eq. (7) can be obtained by minimizing the variational 
functional in Eq. (33) with respect to co^. After setting the first deriv- 
ative of $ in Eq. (33) equal to zero, a system of linear algebraic equations 
is obtained as follows: 


dA CO (34) 


t t 

ACC A 

X— X— 


+ A^C^C 
- 


a ] 


dA = 


. t^t„ 

A C B 


where matrix A is constant for each triangle. Integrating and adding the 

and terms in Eq. (34), a square matrix is obtained for each of the 

triangular elements. The matrix C^^ is also constant for each triangular 
element. Each of the matrices in Eq. (34) is given in Appendix B. 
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Vorticity transport equation - In the solution of the vorticity trans- 
port equation, the variation of velocities over each finite element was ap- 
proximated linearly as 


^ " ‘’1 ^2 ^3 


u. 


u. 


u. 


= B u 


(35) 


^1 ^2 ^3] 


= B V 


(36) 


Substituting Eqs.(32), (35) and (36) into the variational form of the vor- 
ticity transport equation, Eq, (24), the following discrete functional was 
obtained: 


. . 3 w , 

$ = I “ O ^ I 


V to 


B^B + B^B 
— X— X —y—y 


0) 


+ (O^B^B u B w + w^B^B v B to V dA 
X— y — 

I 


(37) 


The functional $ in Eq. (37) was then minimized by differentiating with re- 
spect to w. The resulting system of ordinary differential equations with 
respect to time was written in matrix form as 



3(0 



B^B + B^B 
— X— X — y-^ 


1^1 H. ix 


B^B V B 

yi 


dA to (38) 


The velocity distribution obtained from the solution of Eq. (34) together 
with the instantaneous vorticity distribution is used to calculate the time 
derivative of the vorticity distribution. The coefficient matrix of the 
left-hand side of the equation does not change with respect to time. The 
right-hand side of the equation is calculated at each time step of the nu- 
merical integration. Each of the matrices in Eq. (38) is given in Appendix 
B. 


Equation for pressures - The pressures were also assumed to vary linearly 
over a finite element and were defined as follows: 
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( 39 ) 


Repeating the same analytical procedure, a variational functional was ob- 
tained by substituting Eqs. (35), (36) and (39) into Eq. (17) to obtain: 


$ = 



B^B 

— X— X 


B^B , 
-y-Tj 


E dA - 



fiiEdA 


(40) 


By minimizing the functional in Eq. (40) with respect to nodal pressures, 
the resulting system of algebraic equations can be written as 



+ B^B 


dA £ 



(41) 


The right-hand side of Eq. (41) can be calculated in terms of the var- 
iation of velocities over the finite element from Eq. (11) 


Q 


■ 3x2 3^2 


(42) 


The second derivatives of the stream functions in Eq. (42) can be calcu- 
lated from the definition of ip for a finite element in Eq. (29). By tak- 
ing derivatives of the stream functions in Eq. (29), the second derivatives 
can be written in discrete form as 


Y = J K (43) 

where J is a rectangular matrix as shown in Appendix B. The nodal values 
and the derivatives of stream functions, in terms of area coordinates, can 
be defined in Eqs. (44) and (45) respectively. 


= 


'^xxl *^yyl ^^xyl '^xx2 ^yy2 ^^xy2 ^xx3 '^yy3 ^^xy3 


(44) 





_3jf|_ d^ip 

3 Ci 9 e 2 34953 


'^l'^3 


(45) 
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Substituting Eq. (45) into Eq. (43), the second derivatives of the stream 
functions can be obtained for each element in terms of nodal stream func- 
tions as follows: 


1=11 (46) 

Matrix T is obtained from matrix ^ in Eq. (B6) and matrix ^ A in Eq. (29). 
The resultant matrix T for each finite element is given in Appendix B. 

After the vector ^ in Eq. (41) is defined, the final result of the integral 
on the right-hand side of Eq. (41) was calculated as shown in Appendix B. 


Solution of the System of Equations 


For the numerical solution of the governing equations, each of the Eqs. 
(34), (38) and (41) was written in matrix form including the boundary con- 
ditions. A similar procedure was used for all three equations and the pro- 
gramming was carried out in a uniform manner. 

Stream function solution - The boundary conditions for the solution of 
stream functions in Eq. (34) were included as external constraints in ma- 
trix form as follows: 


c^ = d 
-s — — s 


(47) 


for the free stream boundary and 



(48) 


for the obstacle. Eqs. (47) and (48) can be combined with Eq. (34) to ob- 
tain 


Sc b 


i 


h 

— — s -s 




— s 

c^ 0 0 


6 


d 

— s — — 


— cs 


— s 

b^ 0 0 




0 

— s — — 


bs 




(49) 


In Eq. (49), the c^ and ^ matrices are constant while ^ depends on the 

position and the shape of the obstacle. Separating the time-dependent por- 
tion, Eq. (49) can be written as 
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( 50 ) 


S* b* 


ip* 


h* 

— — s 




— s 

b*^ 0 




0 



— bs 




where 



”s 

c 


'b 


" — 
h 





-s 


— s 


— s 


S* = 



b* = 


It 

4: 


= 



t 


— s 

0 

— s 




c 



d 


6 


— s 

_ 


. . 


— s 


— cs 


(51) 


In this case, is a column vector representing the Lagrange multipliers 

corresponding to the obstacle's boundary conditions. In the case of the 
moving obstacle, matrix ^ is defined at each step of the numerical inte- 
gration. Eq. (50) can be solved to obtain jp* as 

1* = ^ - £5 (^^ £5)“^ (^^ £5) ( 52 ) 

(53) 

! 

(54) 


where 


Hs = 


= S* h" 


r = S* ^ b* 
— s — — s 


matrix is a constant for a particular grid system and has to be inverted 
only once. Different types of obstacles can be defined by specifying the ^ 
vector. Computational efforts are thus reduced considerably. 

Vorticity solution - The boundary condition for the vorticities can be 
written in matrix form as 


c aj = 0 
— w — — 


on the free stream boundary and 


, t • 
b u) 
-w — 


e 

— w 


(55) 


(56) 


on the obstacle boundary. 

The total system of equations for the solution of the vorticity vector 
is then 
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— 


In this case W and matrices are again constant and can be combined as 


where 


The time derivative of the vorticities can be calculated from Eq. (58) as 

w* = q - r (b*^ r )'^ (b*^ q - e ) (6C 

— ^ — w W — W-^ w ^ — w 


where 




= W* ^ h* 
— — w 


r = W* ^ b* 
— w — -w 


Matrix W* is again constant for a particular grid system and matrix ^ is 
defined for a particular obstacle shape. 

Solution for pressures - The boundary conditions for pressure are 


c p = 0 
-T5 ^ - 


at the free stream boundary and 
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( 64 ) 


on the obstacle, 
be written as 


b^ p = 0 

T ^ - 


The total system of equations for obtaining pressures can 


t 

c 


b" 

-V 


c 

-P 


) 

E 


h 


6 

-cp 

= 

0 


[-bpj 


0 


( 65 ) 


Combining matrices £ and c^ 


b* 



E* 


' h*" 


6, 


0 






( 66 ) 


where 


p * = 

1 

1^ 

b* = 

'b ' 
-V 

h* = 

' h ' 
-P 

P* = 

E 


o| 

1 



-IP 

0 

A— 

6 

L-^pJ 


( 67 ) 


The solution of Eq. (66) for p^*is 


E* = 


% 


r (b*^ !■ ) ^ 
-P 


% Hp) 


( 68 ) 


where 


% 


= P*"^ h* 


r = P*”^ b* 
- "T 


( 69 ) 

( 70 ) 


Matrix P^* is assembled and inverted only once for a specified finite element 
gridwork. 
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Detailed Analysis of the Boundary Conditions 


The assembly of the matrices corresponding to the boundary conditions 
can now be described in further detail. The governing equations have to be 
defined in a discrete form suitable for computer applications. 


Stream function analysis - The matrix in Eq. (47) represents the 

boundary conditions applied at the free stream boundary. The stream func- 
tions and the velocities at the free stream are described depending on the 
geometry and flow conditions. Vector ^ represents the obstacle boundary 

conditions as described on page 8. 

Matrix ^ of Eq. (48) has two parts; the first half specifies constant 

stream functions on the obstacle and the second half makes the tangential 
velocity on the obstacle equal to zero. The first half of the matrix con- 
sists of terms for the nodes with equal values of stream functions 




m 


- ijj = 


(71) 


The second half of matrix ^ specifies that the tangential velocities 

are equal to zero again at specific points on the obstacle boundary. The 
tangential velocity u^ was defined from figure 4 in terms of the horizontal 

and vertical velocities as 


u 




- u sin 9 + V cos 0 


(72) 


The angle 0 is defined by the obstacle geometry and can be included into the 
matrix formulation of the boundary conditions as follows: 


u 




C A sin 0 + C A cos 9 
-y- -X- 


1 


(73) 


The discrete coefficients of Eq. (73) were put in matrix form as shown in 
Appendix B. 

Vorticity analysis - The equations specifying zero values of vorticities 
at the free stream boundary were assembled in matrix form. The components 
of the c^ matrix were defined at points on the free stream boundary. As de- 
fined in Eq. (72), the tangential component of the velocity was calculated 
again in terms of the nodal velocities and stream functions. Then matrix b 
contains the area coordinates of each of the points of the obstacle. Col- 
umn vector e is defined at each point by: 


e 

— w 



At 


0) 


(74) 
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where o) is obtained from Eq. (32) and from Eq. (73). Substituting Eq. (6) 

into Eq. (72) and taking derivatives with respect to q, the normal gradient 
of the tangential velocity can be written as 



2 2 

4^ sin 20+4' cos 6 + 4* sin 0 
xy ’^xx yy 


(75) 


The second derivative of stream functions has already been presented in Eq. 
(46). 

Analysis for pressure - Matrix c^ of Eq. (63) , specifying zero pressure 
conditions, is identical to matrix c^ for the vorticity boundary conditions. 
In Eq. (64) matrix ^ is formed to satisfy the boundary conditions on the ob- 
stacle surface by specifying that the normal gradient of the pressure is 
equal to zero. For the angle 0 calculated previously in Eq. (72), the normal 
gradient of pressure is 


9p 3p 

3n " 3y 


cos 0 


- sin 0 
dX 


(76) 


Using the discrete finite element representation of the pressures, Eq. (76) 
can be written as 


3ri 


B cos 0 

-y 


B sin 0 

— X 


E 


(77) 


Each of the three terms in Eq. (77), defining the boundary conditions at one 
point, is inserted into the system of equations corresponding to the nodes of 
the surrounding triangle. 


DISCUSSION AND RESULTS 


The analytical procedure was tested by calculating first the unsteady 
viscous flow past a circular cylinder. The procedure was then applied for 
the analysis of flow around an airfoil. The obtained results and some of the 
difficulties related to the initial and boundary conditions will now be des- 
cribed. Accuracy and stability of the solution and the numerical integration 
technique, along with convergence rates and computer time, will also be 
discussed. 


Initial and Boundary Conditions 


The boundary conditions were employed in a discretized form. Several al- 
ternatives for satisfying the boundary conditions on the obstacle surface 
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were tested using different discretized constraint equations. This investi- 
gation was carried out to decide which set of constraint equations yielded 
the best accuracy. 

For the stream functions, the obstacle boundary conditions were defined 
by specifying that the tangential velocities are equal to zero and that the 
stream functions are equal to each other at several points on the obstacle 
boundary. The same boundary conditions were also defined by specifying that 
the horizontal and vertical velocities are equal to zero at the same points. 
The number and positions of the points at which the boundary conditions are 
specified were tested and compared. It was found that the direct applica- 
tion of zero tangential velocity condition was as accurate as specifying that 
both horizontal and vertical velocities are equal to zero. 

For the analysis of vorticities, the boundary conditions were defined by 
setting the vorticity equal to the normal derivative of the tangential vel- 
ocity at the obstacle boundary. The accuracy of the boundary conditions was 
tested by imposing an additional boundary condition at the leading edge of 
the cylinder and the airfoil at zero angle of attack. When the vorticity 
was artificially set equal to zero at the leading edge, the obtained results 
from the numerical integration were close to the results obtained without 
this additional condition. 

Another important feature of the developed procedure is the description 
of the initial vorticity distribution. At the beginning of the numerical in- 
tegration, the vorticity distribution was assumed to be equal to zero for the 
entire velocity field. With advancing time, vorticities develop on the ob- 
stacle surface due to the vorticity boundary conditions and then diffuse 
throughout tlie flow region. In the initial stages, the disturbances intro- 
duced by these boundary conditions are quite high in magnitude. These dis- 
turbances cause instability in the numerical integration of the vorticity 
transport equation. To eliminate this problem, an underrelaxation factor was 
applied for the first few steps of the numerical integration. A small per- 
centage of the vorticities, generated at the obstacle surface, was included 
in the equations. Once these vorticities diffused away from the obstacle 
boundaries, the underrelaxation factor was removed and the integration was 
continued with the full value of the surface vorticity distribution. 

Another computational difficulty encountered during the solution was re- 
lated to the position of the points representing the obstacle on the finite 
element gridwork. Singularities develop in the solution routine, when the 
nodes describing the obstacle are concentrated in any one finite element tri- 
angle. This problem is a result of the number of boundary conditions spec- 
ified in each finite element. Since the stream functions vary as a third or- 
der polynomial, up to three boundary conditions can be specified in each tri- 
angle. However, since .there are two boundary conditions applied to each ob- 
stacle node in the stream function analysis, only one obstacle node can be 
placed in a single triangle. Similarly, since vorticities and pressures vary 
linearly over each finite element and only one type of boundary condition is 
applied in the solution, again only one node can be considered for each tri- 
angle, If more than one node is specified in a triangle, singularities arise 
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in the matrices shown in Eqs. (51), (59) and (67). A similar problem occurs 
if a boundary condition is duplicated. For example, if the vertical and tan- 
gential velocities are specified as the boundary conditions, these two condi- 
tions may be equivalent at some points on the obstacle. 


Efficiency of the Method of Solution 

The systems of linear algebraic equations, Eqs. (49), (57) and (65), 
were solved by direct inversion of the coefficient matrices shown in Eqs. 
(52), (60) and (68). Since three unknowns are considered at each node of the 
gridwork for the stream function solution, the order of the largest matrix is 
equal to three times the total number of nodes plus the number of boundary 
conditions at the free stream. Although the resulting matrix is relatively 
large, it is well -conditioned and the inverse of this matrix was obtained 
quite accurately by the method of partitioning. The accuracy of the results 
was checked from the symmetric properties of the flow around a circular 
cylinder. 

For a UNIVAC 1108 computer with 65K core storage, the computation time 
needed to invert the largest matrix (458 x 458) was slightly over 12 minutes. 
A significant portion of this time was used up in transferring data between 
the low speed drum and the high speed core. This was necessary because the 
core was able to hold only one-fourth of the total matrix at any one time. 

Since the major portion of the system of equations was solved by the in- 
version of matrices S*, W*, and P*, only the constraint equations represent- 
ing the obstacle's boundary conditions were solved at each step of the numer- 
ical integration. The computer time needed for the solution of several sam- 
ple problems, presented in this analysis, is summarized below: 


Re 

At 

No. Time Steps 

Comp . Time 

40 

1.00 

110 

7 . 5 min . 

1000 

0.10 

120 

8.5 min. 

10000 

0.01 

300 

15.5 rain. 


Again, most of the time used by the computer was for the transfer of data 
from drum storage to core storage. 

An important source of numerical errors in the analysis was due to the 
numerical integration procedure. The accuracy and stability of the integra- 
tion was checked by comparing results from different integration techniques 
and diff eren t integration time steps . Three different -'numerical integration 
routines were tested, and the results were found to be quite accurate. The 
time increment employed for each of these integration methods was propor- 
tional to the inverse of the Reynolds number. At higher Reynolds numbers, 
although it required more steps to reach a steady state, the total number of 
steps corresponded to a shorter real time, due to the convection terms. 
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Figure 5 shows the convergence to a steady state with different time 
increments . 


Discussion of Results 


The results shown in figures 6 through 23 include contour plots of 
stream functions and vorticity and pressure distributions around a station- 
ary circular cylinder and a stationary and oscillating airfoil. 

The sequence of streamline plots in figures 6, 9 and 11 shows the dev- 
elopment of the wake behind a cylinder. As time advances, the wake region 
becomes larger, as does .the angle at which the flow separates. Streamlines 
around a stationary airfoil are shown in figure 15. At the beginning of the 
integration, the flow is attached to the airfoil. As time advances, vortic- 
ities are generated, forcing the flow away from the surface of the airfoil. 
This process continues until the vorticities have moved past the airfoil in- 
to its wake. It should be noted that these results remain symmetric through- 
out the entire integration time sequence, which is one indication of the ac- 
curacy and stability of the solution. At higher Reynolds numbers, unsymmet- 
rical results are also obtained. 

The plots of the vorticity distribution shown in figures 7, 10, 12, 16 
and 18 also show symmetric results. The flow around an airfoil, stationary 
or oscillating, is shown in figures 20 and 22, respectively. Since the ini- 
tial vorticity distribution was assumed to be zero, the first few plots show 
the vorticities located near the obstacle. The following plots show the dif- 
fusion of the vorticities in the surrounding area. 

The pressure contours, or isobars, are shown in figures 8, 13, 17, 19, 21 
and 23. The highest pressure is at the leading edge of the cylinder and the 
airfoil . The surface pressure distribution in figure 14 shows that the pres- 
sure decreases towards the trailing edge until reverse flow occurs. In this 
region, the pressure is of an unsteady nature. 


Conclusions 

The advantages of the developed mathematical procedure will now be sum- 
marized. The analysis of the fluid flow around any arbitrarily shaped obsta- 
cle can be facilitated by specifying a series of points on the boundary of 
the obstacle. The solution routines of the analytical procedure are effi- 
cient. Since the time-dependent boundary conditions were kept separately, 
the remaining portion of the system of equations, which amounts to the bulk 
of computations, was assembled in matrix form and inverted only once. The 
velocity distribution was calculated directly at each time step, without us- 
ing iteration. The major portion of the computational cost involved in the 
analysis was due to core-to-drum transfers of data. Further programming 
efforts to minimize the data handling costs can decrease the computational 
cost considerably. 
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The illustrated results give a good indication o£ the type of flow prob- 
lems for which the method can be applied. The cylinder problem was analyzed 
for testing the method. Numerous solutions for the cylinder already exist. 
The obtained results were compared and good correlation was achieved. A re- 
fined evaluation of the method and the accuracy of the obtained results can 
be achieved through a more rigorous error analysis. The numerical errors 
can be related to the physical problem in terms of the flow characteristics, 
the obstacle shape and the finite element representation of the problem; 
such as the finite element gridwork and numerical integration technique. 

The use of higher order finite elements can also lead to an improvement of 
the present method. 
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APPENDIX A 


AREA COORDINATES 
Definition 

Area coordinates are useful for integration of polynomials over an area 
[12, 13]. In this analysis, area coordinates were employed to approximate 
a variational functional to be integrated over a finite element region. 

Using area coordinates, the position of a point P inside the triangle in 
figure A1 can be specified in terms of the ratio of the area of the subtri- 
angles to the area of the triangle. The ratios for A^^, A 2 , A^ are defined 
as : 



Figure Al. Area coordinate definition for arbitrary point P. 

The C 's are defined as area coordinates. The three area coordinates are 
m 

also related by the expression 


?1 + ?2 ^ 

The area coordinates can be written in terms of rectangular coordi- 
nates by satisfying the conditions at the nodes and using Eq. (A2) as 
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(A3) 





^2 

^ 3 ' 

y 

= 

^1 

^2 

^3 

1 


1 

1 

1 





2 


Inversion of the square matrix in Eq. (A3) leads to the definition of area 
coordinates in terms of the rectangular coordinates: 



■ 


^23 

^32 

^2^3 ^3^2 



II 

^2 

1 

2A 

^31 

^13 

^3yr^i>^3 


y 


^3 


^12 

^21 

Xiy2-X2yi ^ 


1 


where 


X = X - X 

mn m n 


mn 


m 


(A4) 


(A5) 


Differentiation of Area Coordinates 

The derivatives of the area coordinates in x and y directions can be cal- 
culated using Eq. (A4) as 


r 3^1 1 



3x 


^23 


1 


9x 

2A 

>^31 




3x 


^12 


(A6) 
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Integration of the Area Coordinates 

The simplicity of the integration of polynomials in terms of area coor- 
dinates is one of the important advantages of area coordinates. In general 
form the integral can be written as [12]: 



c" rP dA = 2k ^' - P i 

^2 ^3 (m+n+p+2)! 


(A8) 


Eq, (A8] was used extensively in deriving the matrices in Appendix B. 
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APPENDIX B 
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(B6) 

where 


<|) = X y + X y 

1 13 12 213 1 


(j) = X y + X y 

2 32 12 21 23 


(f) = X y + X y 

3 3231 1323 
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APPENDIX C 


FLOW CHARTS 


Several computer programs were prepared in obtaining the presented re- 
sults. The first step of the numerical procedure was to determine the mat- 
rices describing stream functions, vorticities and pressure distributions 
for the finite element mesh without the obstacle; the square matrices , W* 
and jP* shown in Eqs. (51), (59) and (67) were assembled and inverted. Since 
these matrices remain constant, they were inverted only once and then stored 
for use in the solution program. After these matrices were inverted, the 
matrices representing the boundary conditions in Eqs. (52), (60) and (68) 
were added to give the final solution. The overall procedure is summarized 
in figure Cl . 

The matrices defining the finite element mesh were assembled from indi- 
vidual element matrices. These individual matrices were derived from Eqs. 
(B2), (B4) and (BIO) for stream functions, vorticities and pressures re- 
spectively. It should be noted that the order of the element matrix for 
stream functions was three times larger than those for either vorticities 
or pressures, since at each node there were three unknowns; stream function 
and two velocities. For the final matrices S^*, W* and £*, the order of the 
latter two was equal to the number of nodes in the finite element mesh, 
while the order of S^* was three times larger. Since the entire matrix 
could not be stored in the computer at one time, it was assembled one quad- 
rant at a time and stored on magnetic tape as shown in figure C2. The 
smaller matrices, for vorticity and pressure, were assembled at the same 
time in conjunction with the free stream boundary conditions, as shown in 
figures C3 and C4 respectively. 

After the matrices were inverted and stored on tape, the solution was 
obtained using another computer program. This program sets up the matrices 
for the obstacle boundary conditions and the column vectors corresponding 
to the right-hand sides of Eqs. (49), (57) and (65), and then calculates 
the velocity, pressure and vorticity distributions. These distributions can 
be plotted and their development with respect to time can be recorded at any 
time instant. The flow chart for this final program is given in figure C6. 
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Fig. Cl. Overall procedure for the computer solution 












Fig. Cl. Continued. 





Fig. C2. Computational procedure for assembly of the constant matrix for 
stream functions. 
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Fig. C2. Continued. 


2 











Fig. C3. Computational procedures for assembly and inversion of the con- 
stant matrix for vorticities. 
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Fig. C4. Computational procedure for assembly and inversion of the con 
stant matrix for pressures. 











Fig. C5. Inverse routine for the stream function matrix. 
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0^ t j ’ s ^1 


Save this triangle's number 
and these area coordinates 


Continue U 


I 



Do for each 
obstacle node 


Find the triangle for each 
node and its coordinates 


Fig. C6. Flow chart to obtain the final solutions. 
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Fig. C6. Continued. 
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Fig. C6. Continued. 
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Fig. C6. Continued 
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Fig. C7 . Flow chart for the subroutine used in the final program shown in 
fig. C6. 
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Finite element mesh for circular cylinder calculations. 
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Fig. 3. A typical triangular grid element for the finite element analysis 
of the Navier-Stokes equations. 



Fig. 4. Definition of the tangential and normal directions for any point 
on the obstacle. Sines and cosines for point ?2 are defined in terms of 
points and P^. 
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Vorticity, sec.-l Vorticity, sec 



Time, sec. 


Fig. 5a. Stability of the numerical integration of the vor 
ticity transport equation for a point upstream of the ob- 
stacle (Re = 40) . 



Time, sec. 

Fig. 5b. Stability of the numerical integration of the vor 
ticity transport equation for a point near the obstacle 
(Re = 40). 
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circular cylinder (Re = 40). 
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Fig. 7. Vorticity distribution around a circular cylinder (Re = 40). 
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t = 10.5 

Fig. 10b. Vorticity distribution around a circular cylinder (Re = lO’). 
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Fig. 14. Development of the surface pressure distribution 
on a circular cylinder (Re = 40) . 
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Fig. 15. Streamlines around a NACA 0012 airfoil at a 0° angle of attack (Re 




tribution about a NACA 0012 airfoil at a 0° angle of attack (Re 
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Fig. 16b. Vorticity distribution about a NACA 0012 airfoil at a 0° angle of attack (Re = 10 ) 









Fig. 17. Pressure distribution around a NACA 0012 airfoil at a 0° angle of attack (Re 















■essure distribution around a NACA 0012 airfoil 
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21. Pressure distribution around a NACA 0012 airfoil at 
angle of attack (.Re = 10^) 












(t = 0.01 - 


Fig. 23. I’rcssurc distribution about ;»n oscillating NACA 0012 airfoil (kc = 10^“;, 
0.05). (a = 0® - 2.5*). 
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